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ABSTRACT 

We quantify the prospects for using emission lines from rotational transitions of the CO molecule to 
perform an 'intensity mapping' observation at high redshift during the Epoch of Reionization (EoR). 
The aim of CO intensity mapping is to observe the combined CO emission from many unresolved 
galaxies, to measure the spatial fluctuations in this emission, and use this as a tracer of large scale 
structure at very early times in the history of our Universe. This measurement would help determine 
the properties of molecular clouds - the sites of star formation - in the very galaxies that reionize the 
Universe. We further consider the possibility of cross-correlating CO intensity maps with future obser- 
vations of the redshifted 21 cm line. The cross spectrum is less sensitive to foreground contamination 
than the auto power spectra, and can therefore help confirm the high redshift origin of each signal. 
Furthermore, the cross spectrum measurement would help extract key information about the EoR, 
especially regarding the size distribution of ionized regions. We discuss uncertainties in predicting the 
CO signal at high redshift, and discuss strategies for improving these predictions. Under favorable 
assumptions, and feasible specifications for a CO survey mapping the CO(2-l) and CO(l-O) lines, the 
power spectrum of CO emission fluctuations and its cross spectrum with future 21 cm measurements 
from the MWA are detectable at high significance. 

Subject headings: cosmology: theory - intergalactic medium - large scale structure of universe 



1. INTRODUCTION 

In this paper, we consider the possibility of study- 
ing large scale structure at very high redshift (z > 6) 
using rotational emission lines from CO molecules, as 
first proposed by Righi et al. (2008). We further con- 
sider cross-correlating these measurements with upcom- 
ing data from redshifted 21 cm line surveys. In this intro- 
ductory section, we first motivate our study by discussing 
the scientific promise of these measurements. 

One of the primary goals of observational cosmology 
at present is to detect, and elucidate the nature of, the 
Epoch of Reionization (EoR). The EoR is a key stage 
in the history of our Universe when early galaxies and 
quasars turn on and photo-ionize 'bubbles' of neutral hy- 
drogen gas in their surroundings. The ionized bubbles 
grow and merge, eventually filling essentially the entire 
volume of the intergalactic medium (IGM) with ionized 
gas. Detailed measurements of the EoR will help deter- 
mine the properties of the first galaxies and the nature 
of the high redshift IGM, thereby providing important 
observational input for theories of first structure forma- 
tion. 

The most direct, and perhaps the most promising, 
way of studying the EoR is to use the redshifted 21 cm 
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line, which may ultimately provide full three-dimensional 
maps of the reionization process (e.g. Madau et al. 1997, 
Zaldarriaga et al. 2004, Furlanetto et al. 2006a). Moti- 
vated by the promise of this technique, a first generation 
of experiments aimed at detecting the 21 cm signal from 
the EoR is currently underway, including the Murchison 
Widefield Array (MWA) (Lonsdale et al. 2009); the Low 
Frequency Array (LOFAR) (Harker et al. 2010); Preci- 
sion Array for Probing the Epoch of Reionization (PA- 
PER) (Parsons et al. 2010); and the Giant Metrewave 
Radio Telescope (GMRT) (Paciga et al. 2010). These 
first experiments will detect the 21 cm signal only sta- 
tistically, as they will not have the sensitivity to make 
detailed maps of the EoR (McQuinn et al. 2006). The 
surveys will measure, for example, the power spectrum 
of 21 cm fluctuations by binning together many individ- 
ually noisy wave modes. The experiments must over- 
come several observational challenges, especially system- 
atic effects from strong galactic and extra-galactic fore- 
ground emission, coupled with instrumental effects from 
the beam, polarization response, and calibration errors 
(e.g. Liu et al. 2009, Datta et al. 2010, Harker et al. 
2010, Petrovic & Oh 2010). 

Given the statistical nature of the observations, and 
the challenge of the measurements, it is natural to ask 
two related questions. First, how can one ensure that the 
measured signal truly originates from the high redshift 
IGM and is not the result of, for example, residual fore- 
ground contamination? Second, how does one best inter- 
pret the statistical measurement and robustly infer key 
information about reionization? One way of potentially 
addressing both of these questions is to cross-correlate 
the redshifted 21 cm signal with a galaxy survey at very 
high redshift (Furlanetto & Lidz 2007, Wyithe & Loeb 
2007a, Lidz et al. 2009). This measurement is much 
less sensitive to foreground contamination than measure- 
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ments of the 21 cm auto power spectrum since most of 
the anticipated 21 cm foregrounds come from low red- 
shift - primarily galactic synchrotron - and therefore do 
not impact the high redshift galaxy survey. In addition, 
the cross spectrum provides a more direct tracer of bub- 
ble growth during reionization than the auto spectrum: 
the cross spectrum turns over on scales smaller than the 
size of ionized bubbles around the galaxies in the survey 
(Lidz et al. 2009). Hence the cross spectrum can both 
confirm the high redshift origin of a putative redshifted 
21 cm signal, and facilitate its interpretation. 

The difficulty with this proposal is that high redshift 
galaxies are, of course, very difficult to detect. Fur- 
ther, in many respects the galaxy and 21 cm surveys 
are ill- matched (Furlanetto & Lidz 2007): galaxy sur- 
veys typically have high angular resolution but - with- 
out expensive spectroscopic follow-up they have poor 
spectral resolution. On the other hand, the 21 cm sur- 
veys have high spectral resolution, yet poor angular res- 
olution. Further, the purely transverse modes that are 
easiest to measure with a galaxy survey are lost in the 
21 cm data owing to foreground contamination/cleaning. 
Finally, galaxy surveys become increasingly challenging 
at progressively higher redshifts. Nevertheless, it may 
be possible to detect the cross spectrum using upcoming 
surveys for Lyman-alpha emitting galaxies (LAEs) with 
the Hyper-Suprime Cam on the Subaru telescope (Lidz 
et al. 2009). 

In this paper, we consider the complementary possi- 
bility of cross-correlating 'intensity maps' of CO emis- 
sion from high redshift galaxies with upcoming 21 cm 
data. Alongside studies of the 21 cm signal from the 
EoR, it was pointed out that similar experiments could 
map large scale structure after reionization by detecting 
21 cm emission from residual neutral hydrogen locked 
up in galaxies and damped Ly-a systems (Chang et al. 
2008, Wyithe & Loeb 2007b). These works advocate an 
approach dubbed 'intensity mapping', that is rather or- 
thogonal to that of traditional galaxy surveys. In the 
intensity mapping approach, one simply measures the 
large scale variations in 21 cm emission from many in- 
dividually unresolved galaxies across the Universe. This 
lets one measure large scale modes of interest (such as 
the baryon acoustic oscillation scale), without actually 
resolving individual galaxies, which in turn allows for 
an inexpensive experiment. One might wonder if in- 
tensity mapping observations may be usefully performed 
with other (than 21 cm) emission lines, and also at high 
redshifts during the EoR. Indeed, one such possibility 
was considered by Righi et al. (2008), who suggested 
performing an intensity mapping experiment using rota- 
tional emission lines from the CO molecules residing in 
star- forming galaxies at high redshift. 8 This possibility 
is in part motivated by detections of CO emission lines 
from individual quasar host galaxies at z > 6 (Walter 
et al. 2003). A further intriguing possibility is then to 
cross-correlate a high redshift CO intensity map with up- 
coming 21 cm observations. In particular, note that CO 
intensity mapping experiments may be better matched 

8 A few earlier papers suggested performing high redshift inten- 
sity mapping observations with various atomic fine structure lines 
(Suginohara et al. 1999, Basu et al. 2004, Hernandcz-Montcagudo 
et al. 2007, 2008). 



to future 21 cm observations than more traditional high 
redshift galaxy surveys. 

Moreover, a CO intensity mapping experiment also 
provides a unique large-scale view of star formation dur- 
ing the 'cosmic dawn.' Observations from the Hubble 
Space Telescope and ground-based telescopes are pro- 
viding first glimpses into this era, and in the future new 
instruments like ALMA, JWST, and large near-infrared 
ground-based telescopes will provide even more detail. 
However, these instruments (and especially ALMA) are 
all restricted by relatively small fields of view, requiring 
targeted surveys over small scales. Intensity mapping 
complements such efforts beautifully, providing maps on 
several-degree scales. It will therefore provide an unbi- 
ased view of the distribution of CO-emitting gas - and 
hence a window into the distribution of star formation 
over a large range of environments - that is difficult to as- 
semble from targeted measurements of individual galax- 
ies. Together with the detail provided by ALMA, such a 
measurement will provide a complete census of molecular 
emission during the EoR. 

Here we aim to extend the work of Righi et al. (2008) 
and quantify the survey requirements for measuring the 
power spectrum of CO emission fluctuations, as well as 
the 21 cm-CO cross power spectrum. One difference with 
this previous work is that these authors considered the 
angular power spectrum of CO fluctuations func- 
tion of spectral resolution. In this work, we instead treat 
the CO data in a fully three-dimensional way - translat- 
ing observed frequency for gas emitting in a given CO 
line to line of sight distance (modulo peculiar velocities) 
- as planned for future 21 cm data sets. As with the 
21 cm data, one can clean foreground emission from the 
CO data by removing the spectrally smooth component 
of the emission in each spatial pixel. In this paper, we 
focus on two particular CO emission lines: the CO (1-0) 
line and the CO (2-1) line. These two lines are convenient 
because their redshifted emission from the EoR can po- 
tentially be observed from the ground, i.e., the proposed 
observations do not require a costly space-based mission. 
For reference, gas at z = 7 emitting in the CO(1-0) and 
CO(2-l) lines can be observed at v = 14 Ghz and v = 29 
Ghz respectively. Having two CO emission line tracers of 
high redshift structure further guards against foreground 
contamination in the CO data. 

The outline of this paper is as follows. In §2, we de- 
scribe our model for star formation and the specific in- 
tensity of CO emission from high redshift galaxies. §3 
calculates the spatially averaged CO emission from star- 
forming galaxies as a function of redshift. In §4 we then 
use numerical simulations of high redshift structure for- 
mation and reionization (McQuinn et al. 2007a, 2007b) 
to compute the power spectrum of spatial fluctuations in 
the CO emission, while §5 computes the cross spectrum 
between the CO and 21 cm signals. §6 discusses the im- 
pact of foreground emission on the CO measurements, 
and describes the benefit of detecting emission from two 
different rotational lines from emitting gas at a given red- 
shift. §7 then calculates the detectability of the auto and 
cross spectrum signals for various plausible survey config- 
urations. In §8 we compare our results with other related 
calculations, especially closely related work from Carilli 
(2011) and Gong et al. (2011). In addition to these re- 
cent papers and Righi et al. (2008), our work has some 
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overlap with Visbal & Loeb (2010), who recently con- 
sidered cross-correlating various spectral lines (including 
CO lines) along with the redshiftcd 21 cm line. Their 
work focused on studying large scale structure with these 
lines after reionization, and so our work which considers 
the EoR, is nicely complementary. Finally, we conclude 
in §9 and mention possible follow-up work. Through- 
out we consider a ACDM cosmology parametrized by 
n s = l,cr 8 = 0.8, Q m = 0.27, 0\ = 0.73, n b = 0.046, 
and h — 0.7, (all symbols have their usual meanings), 
consistent with the latest WMAP constraints from Ko- 
matsu et al. (2010). 

2. MODELING THE CO BRIGHTNESS TEMPERATURE 

In order to model the CO emission from high redshift 
galaxies, we adopt a simple model to connect the strength 
of this emission to the abundance of the dark mater halos 
that host CO luminous galaxies. In brief, we first relate 
a galaxy's star formation rate (SFR) to the dark matter 
mass of the galaxy's host halo. We then assume that 
empirical correlations between a galaxy's SFR and its 
luminosity in the CO lines of interest, measured from 
galaxies at z < 3, apply also at higher redshifts during 
the EoR. Throughout, we consider how uncertainties in 
this model impact our theoretical predictions. 

First, let us consider the specific intensity of the CO 
emission. If the proper volume emissivity of CO emitting 
galaxies at redshift z', emitting at an observed frequency 
of fobs, is e [^obs(l + z')], then the specific intensity of 
the observed CO emission can be determined by solv- 
ing the (absorption free) cosmological radiative transfer 
equation. The solution at z = and observed frequency 
t'obs is: 



I(lSobs,Z = 0) 
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,0 dz' (1 
In this equation dl/dz' denotes the proper line element 
dl/dz = c/[{l + z)H(z)}. 

Next we would like to compute the proper volume emis- 
sivity of the CO emitting sources. For simplicity, we as- 
sume that a (halo mass independent) fraction, /duty, of 
halos with mass larger than M COjm j n actively emit in the 
CO line of interest at any given time. Within the ha- 
los that host CO luminous galaxies at a given time, we 
further assume that the CO luminosity is directly pro- 
portional to the halo mass. Finally, we approximate the 
profile of each CO emission line as a delta function in 
frequency. The specific luminosity is then given by 

L v = A5d{v ~ vj)M, (2) 
where Sd denotes a Dirac delta function, vj is the rest 
frame frequency of the transition of interest, M is the 
mass of the dark matter host halo, and A is a propor- 
tionality constant that we will discuss shortly. Denot- 
ing the (co-moving) halo mass function at redshift z by 
dn co (z)/dM, the proper volume emissivity is: 

e(v, z) = A5 D {v - vj){1 + z) 3 /dut y x 
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Inserting Equations 2 and 3 into Equation 1 we find: 
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In this equation is \, s = uj/(1 + zj) is the observed fre- 
quency for gas at redshift zj emitting in a CO line with 
rest frame frequency vj. The equation implies that, in 
our simple model, the specific intensity of CO emission 
is directly proportional to the fraction of matter in halos 
of mass greater than M CO:m ; n . We will generally quote 
results for the CO emission in terms of an equivalent 
brightness temperature, computed using the Rayleigh- 
Jeans Law, T co Kbs) = c 2 I(v ohs )/(2k B vl hs ). 

2.1. The Star Formation Model 

Next we need to specify the parameter A, which relates 
CO luminosity and host halo mass. The first ingredient 
involved in specifying A is to connect the star formation 
rate and host halo mass, M. Here we adopt the simplest 
plausible model and assume that the baryon fraction in 
the host halo follows the universal value, Q(,/f2 TO , and 
that a fraction /* of these baryons are turned into stars 
at a constant rate over a time scale t s . We assume that 
/* and t s arc themselves independent of host halo mass. 
In this case the SFR is (Loeb et al. 2005): 

SFR = h^^f = O.^Moyr- 1 x 



0.1 









10 8 yrs 






0.17 




t s 





M 



l0 9 Mr. 



(5) 



We assume /* = 0.1 and t s — 10 8 yrs. These values are 
consistent with those found by Stark et al. (2007) to 
match the luminosity function of Lyman Break Galaxies 
(LBGs) at z = 6 in a very similar model. Similarly, by 
comparing the UV luminosity function and the 2 point 
correlation function, Lee et al (2009) find the typical du- 
ration of star formation is < 0.4Gyr at z w 4—6. The star 
formation time scale t s sets the duty cycle for star for- 
mation activity to be t s /tji', here tn is the Hubble time 
at the redshift of interest which is also approximately 
the 'age' of the dark matter halos. The duty cycle for 
CO luminous activity, /duty, may differ from t s jtu but 
is likely comparable to this if CO is excited by starbursts 
(see §2.2). Our fiducial model assumes /duty = 0.1, which 
is comparable to t s jtn at the redshifts of interest. 

It is also useful to consider the star-formation rate den- 
sity (SFRD), p\(z), in our model. This depends on the 
collapse fraction of halos with mass larger than M s f >m ; n , 
where M s f im i n denotes the minimum dark matter halo 
mass that hosts a star-forming galaxy. Note that M s { im i n 
may be smaller than M co m i n if some galaxies form stars 
yet do not emit in the CO line (see §2.2). The SFRD is 
then given by: 



p*(z) = O.O79M yr" 1 Mpc~ 



0.1 



t H (z = 7) 



t H (z) 



/coll ( Msf ,min, z) 



0.1 



.(6) 



(4) 



The model star formation rate in the above equation de- 
pends on the Hubble time at the redshift of interest; the 
number above is scaled to the Hubble time at z — 7. Here 
/coii(A4f,min, z) denotes the fraction of mass in halos of 
mass greater than M s f imin at redshift z; it is the halo 'col- 
lapse fraction'. The value of 10% for the collapse fraction 
in the above equation gives its order of magnitude at the 
typical redshifts and halo masses of interest. 
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Following Carilli (2011), it is useful to compare the 
SFRD in our model to the critical star formation rate 
density required to overcome recombinations, and keep 
the Universe ionized (Madau et al. 1999). After reion- 
ization completes, the star formation rate density must 
clearly exceed this critical value, while during most of 
the EoR one would expect the SFRD to be comparable 
to this critical rate. Hence although our star formation 
model is simplistic, and the parameters in Equation 5 
are uncertain, it is nevertheless constrained by requiring 
the star formation density to be comparable to the crit- 
ical value. If pi,(z) were much less than (greater than) 
the critical rate, the Universe would be drastically under- 
ionized (over- ionized) at the redshift of interest. 

This critical SFRD depends on three uncertain quanti- 
ties: the number of ionizing photons produced per baryon 
converted into stars (which in turn depends on the stel- 
lar IMF and the spectral shapes of the stars in the ion- 
izing galaxies), the fraction of ionizing photons that es- 
cape from the host halo and ionize the IGM, and on the 
dumpiness of the IGM. Adopting the ionizing spectrum 
expected for a Salpeter IMF, solar metallicity, and plau- 
sible values of the escape fraction and clumping factor, 
the critical rate is given by (Madau et al. 1999, Munoz 
& Loeb 2010): 



P*,crit = O.OSOMoyr^Mpc- 
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The above choice of clumping factor, C — 5, is similar 
to that measured from recent hydrodynamical simula- 
tions (e.g. Pawlik et al. 2009, McQuinn et al. 2011), 
although there are still theoretical uncertainties in this 
estimate. In particular, note that the clumping factor 
that enters Equation 7 is that of the ionized gas in the 
IGM. The clumping factor then depends on which gas is 
considered to be part of the IGM, and which gas belongs 
to a galaxy, as well as the properties of reionization it- 
self, such as the spatial distribution and spectrum of the 
ionizing sources. The value of the escape fraction in the 
above equation, / csc = 0.1, is similar to that measured 
from LBGs at z <~ 3 (Nestor et al. 2011). However, the 
structure of the faint dwarf galaxies that likely reionized 
the Universe may be rather different than that of these 
LBGs, and the reionizing galaxies may hence have rather 
different escape fractions than the z <~ 3 LBGs. While 
the uncertainties are substantial, Equation 7 still pro- 
vides a useful indication of plausible values for the high 
redshift star formation rate. 

The claim that the SFRD should not differ greatly from 
this critical value at the redshifts of interest is strength- 
ened by considering measurements of the mean trans- 
mitted flux in the Ly-a forest at high redshift. These 
suggest that the total emissivity of ionizing photons at 
z = 5 — 6 is comparable to the critical value required to 
keep the Universe ionized (Miralda-Escude 2003, Bolton 
& Haehnelt 2005, Faucher-Giguere et al. 2008, McQuinn 
et al. 2011). In addition to the Ly-a forest measure- 
ments at z ~ 5 — 6, WMAP constraints on the electron 
scattering optical depth (Komatsu et al. 2010) suggest 
significant amounts of star formation at z > 6; it is hence 
unlikely that the SFRD falls steeply right above z~5 — 6. 
The Ly-a forest data indicate a hydrogen photoioniza- 
tion rate of Thi « 5 x 10 _13 s _1 at z = 5.5 (see the 
above papers for a discussion of the uncertainties in this 



measurement). Combining this with a recent constraint 
on the mean free path to HI ionizing photons, Ahi, from 
Songaila & Cowie (2010), and assuming that the ionizing 
spectrum is a v~ 2 power law near the HI photoionization 
edge gives: 

Thi 

P*Xy- 



= O.O39M yr _1 Mpc~ 



5 x lO-^s" 1 
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In this equation, i\T p h is the number of ionizing pho- 
tons emitted per second for an SFR of 1M Q yr -1 , and 
10 53 1 s _1 is the value expected for a Salpeter IMF and 
solar metallicity (Madau et al. 1999). 9 The num- 
ber in brackets for Ahi is the preferred value of the 
mean free path in units of proper Mpc from Songaila 
& Cowie (2010). Taken together, Equations 7 and 8 
argue that p*(z) is unlikely to differ significantly from 
P*,crit for most of the redshifts of interest, except per- 
haps in the earliest phases of the EoR, at which point 
the CO measurement is likely impossible anyway. Fur- 
ther, the equations quantify the expected SFRD, which 
is p*(z = 7) ~ O.O5M yr -1 Mpc -3 , although with sig- 
nificant uncertainties from the unknown value of C/ f esc . 

Finally, we compare our model SFRD with these num- 
bers. The model depends on M s f im ; n , the minimum mass 
halo hosting a star-forming galaxy. One plausible value 
for M s f m i n is the halo mass corresponding to a virial 
temperature of T vir = 10 4 K - M sfimin w 10 8 M Q at 
the redshifts of interest - above which gas is able to 
cool by atomic hydrogen line emission, condense, and 
form stars (e.g. Barkana & Loeb 2001). This mass 
scale is sometimes called the 'atomic cooling mass'. For 
M s f im i n = 1O 8 M , we compute the model SFRD of Equa- 
tion 6 at z = 7 using the Sheth-Tormen (2001) halo mass 
function, finding p*(z = 7) = 0.060M o yr" 1 Mpc~ 3 . 
This is indeed comparable to the critical rate (Equations 
7, 8). Feedback from photoionization and supernovae ex- 
plosions may, however, limit the efficiency of star forma- 
tion in these small mass halos, raising M s f jm i n , although 
the precise impact of these effects is still unclear (e.g. 
Dijkstra et al. 2004). If M sf . min = 1O 9 M , for exam- 
ple, the model SFRD drops to p+(z = 7) = O.O29M 
yr -1 Mpc -3 , just a little lower than the critical rate for 
'/csc = 0.1, C = 5 (Equation 7). The SFRD in this 
model is closer to the critical rate if the escape fraction 
is higher, the clumping factor lower, or if the normaliza- 
tion of the SFR-M relation is higher (Equation 5). With 
the latter solution, the SFRD would again be similar to 
that in our atomic cooling mass model. In other words, 
normalizing the model to the critical SFRD removes at 
least some of the freedom from varying M s f im ; n (see also 
Carilli 2011, and the present paper, §4.3). Note also that 
-Msf,min might be smaller than the atomic cooling mass 
if molecular hydrogen cooling is efficient despite negative 
feedback from dissociating UV radiation (Haiman et al. 
1997). Given these and other uncertainties, we will sim- 
ply use the atomic cooling mass as our fiducial M s f :m ; n , 
along with the SFR-M relation of Equation 5, since this 
model produces approximately the expected SFRD. 

9 For simplicity of presentation, we do not include separate 
bracketed terms illustrating an additional (weak) dependence on 
the assumed shape of the ionizing spectrum. 
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2.2. CO Luminosity 

The next crucial, but least certain, ingredient in our 
model is to relate the CO luminosity and star formation 
rate at high redshift. We will consider this from both 
an empirical point of view and from a theoretical per- 
spective. Observationally, CO emission lines have been 
detected all the way out to z ~ 6 from the host galaxies 
of bright quasars, providing part of the impetus for our 
present study (Walter et al. 2003, Wang et al. 2010). 
Similarly, CO emission lines have been detected at high 
redshift (z ~ 4.5) from submilimeter galaxies with ex- 
treme starbursts (SFR ~ 1O 3 M yr _1 ) (e.g. Schinnerer 
ct al. 2008). Much less clear observationally is whether 
more normal galaxy populations, without extreme star- 
bursts, are CO luminous at high redshift: so far there 
are only upper limits on CO(1-0) emission from a pair 
of z = 6.6 LAEs (Wagg et al. 2009). From a theoret- 
ical perspective, the CO emission depends on complex 
astrophysical processes that are difficult to model from 
first principles. For example, it depends on the struc- 
ture of high redshift galaxies and the spatial distribution 
of their star forming gas, and the metallicity and tem- 
perature of the molecular clouds in these galaxies. The 
temperature of the molecular clouds in turn depends on 
the level of heating from the CMB, photoelectric heating 
by dust grains and heating from starbursts, AGN activ- 
ity, supernova shocks, and cosmic rays, as well as cooling 
by atomic fine structure lines and molecular lines (e.g. 
Obreschkow et al. 2009). 

Indeed, Obreschkow et al. (2009) present a detailed 
model for CO emission from high redshift galaxies, in- 
corporating models for many of these physical processes. 
Their model starts from a semi-analytic model of galaxy 
formation placed on top of the Millenium simulation 
(Springel et al. 2005), and then incorporates additional 
modeling to split the cold gas into atomic and molecular 
components and for converting H2 abundances into CO 
luminosity. One approach would to use the Obreschkow 
et al. (2009) model directly, as done in Gong et al. 
(2011). While we draw heavily on the discussion in 
Obreschkow et al. (2009) to inform our modeling, we 
prefer a different approach for several reasons. First, the 
minimum resolved halo mass in the Millenium simula- 
tion is 3 x 10 10 M Q ( Springel et al. 2005), more than 
two orders of magnitude larger than the atomic cool- 
ing mass; hence the hosts of the likely ionizing sources 
arc not resolved in this calculation. Furthermore, the 
models for these massive galaxies make assumptions: 
e.g., that all galaxies are quiescent, virialized exponen- 
tial disks — which drives the partition between H2 and 
HI — which is unlikely to apply in the small, constantly 
merging systems which drive the CO signal at high red- 
shift. Next, observed galaxy properties at low and in- 
termediate redshift are relatively well understood — and 
the semi-analytic models are well-calibrated at these red- 
shifts - but the properties of low luminosity galaxies at 
z > 6 are much less certain. The use of a highly com- 
plex, multi-parameter model with potentially inapplica- 
ble physics obfuscates the key requirements-and associ- 
ated uncertainties — for a detectable signal. As a result, 
we will adopt a much simpler model for CO luminosity 
and try to identify key sources of uncertainty. 

Let us start by considering empirical correlations be- 



tween star formation rate, far-infrared luminosity, and 
CO luminosity as measured in relatively nearby galaxies 
(see also Carilli 2011). Wang et al. (2010) fit a correla- 
tion between far-infrared luminosity and CO luminosity 
using measurements from galaxies at z — — 3, includ- 
ing spiral galaxies, LIRGs, ULIRGs, and SMGs. Specif- 
ically they fit a power law relation between far-infrared 
luminosity and the velocity integrated CO(1-0) bright- 
ness temperature multiplied by the source area. The lat- 
ter quantity is denoted by ^co(i-o)' an( ^ ^ as dimensions 
of K km/s pc 2 . It is related to the total luminosity of 
the emission line, £00(1-0)7 by the relation: 



L 



CO(l-O) 



1.04 x lO^Lp 



© 



L' 



CO(1-0) 



3.25 x lO^kms-ipc 2 



The Wang et al. correlation is: 



log 



/ Lfir \ 



V L 



1.671og 



CO(l-O) 







Kkms 1 pc 2 



(9) 



4.87, (10) 



with Xfir denoting the far-infrared luminosity. This 
can then be combined with the Kennicutt (1998) cor- 
relation between SFR and Lfir, which is SFR = 
1.5 x lO~ lo M Q yr -1 (L F iR/L0). Assuming this relation 
between far-infrared luminosity and SFR, the Wang et 
al. (2010) fit covers a wide range of SFRs from ~ 
0.1 — 1,000-Mq yr _1 . Combining the Kennicutt relation 
with Equations 9 and 10, we arrive at a relationship be- 
tween the luminosity in the CO(1-0) line, and the star 
formation rate: 

3/5 



ico(i-o) = 3-2 x 10 4 L 



SFR 



Note that according to this relation the CO luminosity is 
a sub-linear function of the SFR. This acts to weight the 
contribution of galaxies with low star formation rates to 
the emissivity (Equation 3) most heavily. This implies 
that CO emission may be especially bright if the scaling 
continues to apply at high redshifts (z > 6) and low SFRs 
(SFR < 0.1M© yr" 1 ), since abundant dwarf galaxies 
with low SFRs are likely ionizing sources. 

One possible theoretical explanation for the sub-linear 
scaling of CO luminosity with SFR was proposed by 
Narayanan et al. (2010): the CO(1-0) transition has 
a low critical density (to thermalize) and so emission in 
this line may be independent of SFR and simply pro- 
portional to the gas density, Lco(i-o) P- Combining 
this with the Schmidt (1959) law, SFR oc p 3 / 2 , gives 
approximately the observed sub-linear scaling. If this 
explanation is correct, it may hold at higher redshift as 
well although as discussed below, a sufficiently strong 
starburst is required to produce a high enough excita- 
tion temperature, T ox , for CO to be observable against 
the high redshift CMB. Lines from significantly higher 
rotational states will have higher critical densities, and 
the excitation rate of these lines will then depend on the 
SFR, and give a different, strong ding of Lqo w / 
SFR, as found observationally (Naryanan et al. 2010). 
These authors' simulations give, however, a similar sub- 
linear scaling for CO(2-l) and CO(1-0), and so this scal- 
ing appears reasonable for the two emission lines of in- 
terest for our study. 
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Let us now consider the various effects that may mod- 
ify the empirical CO luminosity-SFR relation of Equation 
11 at high redshift. One significant worry is that the CO 
emission must be detected against the CMB, whose tem- 
perature grows with redshift as Tcmb oc (1 + z). In the 
limit that the CO excitation temperature, T ex , is equal 
to the CMB temperature, the CO emission is completely 
undetectable at the redshift of interest (Combes et al. 
1999, Obreschkow et al. 2009). As mentioned in the 
beginning of this section, a number of processes may be 
important in determining the CO excitation tempera- 
ture, and it is unclear how these will impact the molec- 
ular gas in the high redshift dwarf galaxies of interest. 
In particular, although high redshift dwarfs are small in 
size and have low luminosities, the relevant quantity for 
computing their excitation temperature is the star for- 
mation density. Specifically, in the optically thick limit, 
T 4 x w Tq MB + T^ oat , where Th ea t reflects the dust heat- 
ing due to starburst and/or AGN activity; in order for 
the CMB not to significantly dilute the CO emission, 
Theat > Tcmb- The molecular gas is likely very dense 
in these high redshift galaxies, and so their star forma- 
tion densities should be large, provided their star forma- 
tion obeys a Schmidt-type law. Note that in starbursts, 
T cx > 40K in observed CO emitters (Obreschkow et al. 
2009), which still gives plenty of contrast with the CMB 
at z ~ 8, for example. Indeed, given the potentially 
higher densities in high redshift galaxies, it is possible 
that high redshift dwarf galaxies are more luminous in 
CO than expected from the low redshift empirical con- 
straints (although note that this is more likely to boost 
detectability in high J transitions). 

A second possible worry is that the metallicity will de- 
crease towards high redshift and that this will rapidly 
diminish the CO emissivity at early cosmic times. This 
is likely less severe a problem than it might first ap- 
pear. First, while the spatial average metallicity is likely 
quite low at z > 6, rapidly star- forming regions will typ- 
ically already be enriched by several generations of stars. 
Second, the low order CO transitions likely come from 
optically thick regions, in which case the emission may 
be nearly independent of metallicity. However, there is 
likely less dust in these high redshift star-forming re- 
gions than in their low redshift counterparts, and the 
dust helps shield CO from UV dissociating radiation (e.g. 
Obreschkow et al. 2009). 

Given the uncertainties, we will stick to the sim- 
plest plausible model outlined in §2. This model is 
parametrized by Af co . m ; n , the minimum host halo mass 
for CO luminous halos; the luminosity of such halos, 
Leo (-^co, min) (with Leo oc M for larger mass halos); 
and the duty cycle, /duty, for CO activity. The latter pa- 
rameter is entirely degenerate with Lco(-^Lo,min) as far 
as the mean specific intensity of the CO emission is con- 
cerned (Equation 4), but the duty cycle further impacts 
the level of Poisson fluctuations in the CO emission (§4). 
The parameter M COim ; n will crudely mimic the scenarios 
discussed above in which low SFR galaxies are dim in 
the CO lines of interest at high redshift. For instance, 
Aico,min niay be large if dwarf galaxies with low SFRs do 
not boost the gas temperature and the excitation tem- 
perature of the low-J CO lines, T ex , enough to be visible 
above the CMB background at z > 6, or if the metal- 
licity of the star-forming regions in these galaxies is too 



low and CO is dissociated result. 

Our aim is to understand which values of these pa- 
rameters produce a detectable CO power spectrum and 
21 cm-CO cross power spectrum signal. In our fidu- 
cial calculations, we assume M CO!m j n = M coo i ~ 1O 8 M , 
/duty = 0.1. We further adopt the (lower redshift) em- 
pirical CO luminosity-SFR relation of Equation 11. For 
simplicity, we take Leo oc M cx SFR and so our model 
will not precisely match the sub-linear scaling of Equa- 
tion 11 for all model galaxies. However, we normalize 
the model CO luminosity to match that demanded by 
the sub-linear scaling of Equation 11 for star forming 
halos with the atomic cooling mass. In other words, 
we fix Lqo(M = Meed) using Equation 11 and assume 
Leo oc M for larger mass halos. In §4.3 we check the 
impact of relaxing this assumption, and generalizing the 
dependence to Leo oc M a . There we find that assuming 
the linear scaling for masses above M coo i, rather than 
the sub-linear one, has little impact on our results. This 
occurs because the emissivity in these models is domi- 
nated by the sources in the smallest halos. For instance, 
the mean emissivity is only a factor of two higher in the 
linear case. 

Combining Equation 5 and Equation 11, we obtain 
Leo = 2.8 x 1O 3 L for halos of mass M = 10 8 M©. Our 
fiducial model then assumes 



Lco(M) = 2.8 x 10 3 L Q 



M 



10 8 M n 



(12) 



The empirical relation (Equation 11) is for the CO (1-0) 
line, and so the above equation strictly applies only to it, 
and not to emission lines from higher rotational levels. 
However, given the uncertainties in the normalization of 
the CO luminosity at high redshift, we will assume that 
it applies to the CO (2-1) line as well, which is very likely 
a conservative assumption. For example, in the limit 
that both lines are optically thick and fee^ex >> hi/j, 
Lco(2-i) = 8Lco(i-o) ( e -g- Obreschkow et al. 2009), 
although the ratio of the luminosities in the two lines is 
unlikely as large as in this limit. 10 The relation in Equa- 
tion 12 then fixes our fiducial choice of the parameter A 
in Equations 2, 3, and 4. 

For the CO(2-l) line, the normalization of Equation 12 
is a factor of 6 higher than that adopted in Righi et al. 
(2008) at the star formation rate of atomic cooling mass 
halos in our model. Our CO(1-0) luminosity is higher 
than theirs by a still larger factor. Their lower normal- 
ization is driven by their matching to M82, which has an 
SFR of <~ 10M Q yr _1 . This is significantly larger than 
the expected SFRs of the low luminosity dwarf galax- 
ies that likely reionized the Universe, and so we prefer 
our normalization which may provide a better match to 
the low luminosity galaxies of interest. In §8 we provide 
further comparison with previous work. 

3. SPATIALLY-AVERAGED CO EMISSION 

We now have all of the basic ingredients necessary 
to calculate CO brightness temperatures in our model. 
Combining Equations 4, 5, and 12 with the Rayleigh- 
Jeans relation between brightness temperature and spe- 
cific intensity, and inserting our fiducial model parame- 

10 For reference, huj jkg = J5.5K. 
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ters, we obtain: 

Tco (v bs) = 2.1/j.K 



/*coll(^^co,min5 %j) 



0.1 



/duty 



o.i 



"2" 


3 


~1 + ZJ~ 


1/2 




J 




8 




1 



(13) 

In this equation, gas emitting at a redshift zj in a CO 
transition of rest frame emission frequency vj is observed 
at a frequency f bs = + Zj). The frequency emit- 

ted in a transition between rotational levels J — > J — 1 is 
j/ j = JvcOi with z^co = 115 Ghz. The fiducial value for 
J in the above equation is J — 2, i.e. the numbers are 
for the CO(2-l) transition. CO gas at the same redshift 
emitting in the CO (1-0) transition will be observed at 
half of the frequency of the CO(2-l) line. This equation 
makes the conservative assumption that the CO (2-1) line 
and the CO (1-0) line have the same luminosity-halo mass 
relation. In the high temperature and optically thick 
limit mentioned in the previous section, and adopting 
Equation 12 for the CO(l-O) line, the brightness temper- 
ature in the CO(2-l) line would be a factor of 8 higher 
than implied by Equation 13. For simplicity, we will 
focus most of the remaining discussion on the CO(2-l) 
line and adopt Equation 13: the reader can rescale the 
results for the CO(l-O) line based on Equation 13 and 
their preferred luminosity-halo mass relation. This com- 
pact equation then allows us to predict the CO bright- 
ness temperature for a region in which a fraction / co u of 
matter has collapsed into dark matter halos hosting CO 
luminous galaxies. Note that Tco and / co ii will fluctuate 
from place-to-place across the Universe, but we suppress 
the spatial coordinate in our notation here. 

We start by computing the volume-averaged CO 
brightness temperature. In order to illustrate some of 
the dependence on model uncertainties, we vary the min- 
imum halo mass hosting a CO luminous galaxy widely, 
between M cc , min = 10 8 , 10 9 , 10 10 , and 1O U M . In 
our model, these halos host star formation rates of ~ 
0.02,0.2,2, and 20M© yr _1 , respectively. Note that this 
minimum mass may be distinct from the minimum mass 
dark matter halo in which stars form (we denote this lat- 
ter mass by M s f im i n ) if some galaxies are simply not CO 
luminous, as discussed in the previous section. 

The results of the mean brightness temperature calcu- 
lation are shown in Figure 1. For this calculation, we use 
the Sheth-Tormen (2001) mass function to calculate the 
collapse fraction, since this formula matches the abun- 
dance of simulated dark matter halos at high redshifts 
fairly well (e.g. Zahn et al. 2007). In all models, the 
mean CO brightness temperature is a steep function of 
redshift, which results because the host halos are rare 
objects at the initial redshifts of interest, and the col- 
lapse fraction consequently grows rapidly towards lower 
redshifts. This steep redshift evolution is enhanced for 
the large values of M CO:mm , in which case the CO lu- 
minous host halos are correspondingly rare. At z ~ 7 
the mean brightness temperature is of order IfiK for 
M co>m i n = 10 8 Mq, but the results fall with increasing 
-^co,min- For example, the mean brightness tempera- 
ture at z ~ 7 is a factor of w 2 smaller than this if 
-Mco.min = 10 9 Mq, and it is a factor of ~ 7 smaller than 
this if M COiroin is as large as 1O 1O M (SFR = 2M© yr" 1 in 
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Fig. 1. — The global mean CO brightness temperature as a func- 
tion of redshift. The curves show the volume- averaged CO bright- 
ness temperature model of Equation 13 for several different values 
of M, 



co, mm i 



the minimum halo mass hosting a CO luminous galaxy. 
The models assume that the minimum host halo mass of star form- 
ing galaxies is the atomic cooling mass, M a { m ; n fa 1O 8 M0, and so 
the curves with larger M co m ; n describe models in which galaxies 
with low SFRs are not CO luminous. If we had instead varied 
A/ s , 



while fixing both M COjll 



M B f m ; n and the SFRD at a 



given redshift, the model variations would be considerably smaller 
(see text). 



our model). The mean emission drops off still further in 
the extreme case that M co mm = lO 11 A/ , in which case 
(Tco) is about a factor of 50 smaller than in the atomic 
cooling mass model. The strong dependence on Af C o,min 
arises because the collapse fraction depends strongly on 
the minimum host halo mass for these rare halos. At still 
higher redshifts, the relevant dark matter halos are rarer, 
and the dependence on M COim j n is still stronger. As we 
will detail subsequently, the dependence shown here has 
potentially strong implications for the feasibility of CO 
intensity mapping experiments at high redshift. 

It is important to emphasize, however, that the most 
important source of uncertainty is not so much in the 
minimum mass of halos that host star-forming galaxies 
(-^sf,min), but rather in the minimum host halo mass of 
CO luminous galaxies (M co m i n ). This is because, re- 
gardless of M s f jm i n , the SFRD must be comparable to 
the critical rate (Carilli 2011) once reionization is signif- 
icantly underway. In the context of our model, boosting 
M s f jln j n necessitates increasing the SFR-M normalization 
to ensure that the SFRD is still comparable to the critical 
rate. In our fiducial model, both (Tco) and the SFRD 
are proportional to the collapse fraction, with (Tco) scal- 
ing with the SFR-M normalization to the 3/5th power 
(because Leo scales sub- linearly with SFR). Hence nor- 
malizing to a fixed SFRD at a given redshift will remove 
most of the dependence on M s f inun at that redshift. We 
quantify this further in §4.3. The only caveat with this 
argument is that the precise normalization of the SFRD 
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is subject to considerable uncertainties (§2.1). Nonethe- 
less, the main concern is not in the value of M s f im ; n but 
that the galaxies that produce most of the ionizing pho- 
tons - likely abundant sources with low star formation 
rates - may be CO dim for the reasons discussed in §2.2. 
The value of M s f im i n does of course impact the strength 
of brightness temperature fluctuations, as we discuss in 
§4- 

Finally, let us discuss how the spatial mean brightness 
temperature depends on the duty cycle for CO luminous 
activity, /duty First note that the SFRD is indepen- 
dent of the starburst duty cycle: as the star formation 
timescalc increases, an increasing fraction of gas needs to 
be converted into stars to produce a given SFRD (Equa- 
tions 5-6). Because of this, the duty cycle impacts (Too) 
most strongly if the duty cycle for CO luminous activity 
is significantly larger or smaller than that for star forma- 
tion. Since CO lines are excited by starburst activity, we 
generally expect the two duty cycles to be comparable, 
as assumed in our fiducial model. If the CO duty cy- 
cle is instead larger (smaller) than this fiducial starburst 
value, the mean brightness temperature will be increased 
(decreased) in direct proportion. 

4. CO BRIGHTNESS TEMPERATURE FLUCTUATIONS 

Unfortunately, it is likely impossible to observe the 
global mean CO brightness temperature directly. The 
difficulty with observing this signal is that there will be 
other sources of emission at the relevant frequencies that 
arc much brighter than the rcdshiftcd CO signal. How- 
ever, similar to the case of the redshiftcd 21 cm signal, 
the CO signal should have structure in frequency space 
while foreground sources will be spectrally smooth (see 
§6). One can use this fact to isolate spatial fluctuations 
in the CO brightness temperature, and detect them ro- 
bustly. 

4.1. Numerical Simulations 

In order to model spatial fluctuations in the CO emis- 
sion and the redshifted 21 cm signal, we use outputs from 
the numerical simulations of reionization of McQuinn 
et al. (2007a, 2007b). These simulations treat the ra- 
diative transfer of ionizing photons in a post-processing 
step performed on top of an evolved N-body simulation. 
The N-body simulation was run using an enhanced ver- 
sion of Gaciget-2 (Springel 2005), and follows 1024 3 dark 
matter particles in a box with a co-moving side length 
of 130 /i" 1 Mpc. The simulation directly resolves halos 
down to ~ 1O 1O M0, but smaller mass halos down to the 
atomic cooling mass scale of M coo i ~ 1O 8 M are incorpo- 
rated with the appropriate abundance and clustering as 
in McQuinn ct al. (2007a). This is done using a merger 
tree algorithm similar to the PThalo code (Scoccimarro 
& Sheth 2002), a widely used code for rapidly generating 
mock galaxy surveys. In the previous section we found 
that halos close to the atomic cooling mass, if indeed CO 
luminous, may provide the dominant contribution to the 
spatially-averaged CO brightness temperature and so it 
is important to accurately capture the impact of these 
halos. Note that for the CO calculations we only require 
simulated halo catalogs from the N-body simulation, and 
not the full reionization simulations. To calculate the 21 
cm field, and the 21 cm-CO cross power spectrum, we use 
reionization calculations from the same fiducial model as 



in Lidz ct al. (2008). In this model, a source's ionizing 
luminosity is directly proportional to its host halo mass 
for halos above the atomic cooling mass. 

4.2. The CO Auto Power Spectrum 

In order to characterize the spatial fluctuations in 
the CO emission, we focus on the power spectrum of 
CO brightness temperature fluctuations. We construct 
(three-dimensional) maps of CO emission from simulated 
halo catalogs using Equation 13 for various redshift out- 
puts and values of M co<m i n . For any given redshift out- 
put, we randomly select a fraction /duty of the simulated 
halos to be CO-luminous, and assume that the remain- 
ing fraction, 1 — /duty, of the halo catalog are not ac- 
tively emitting in CO. We then measure the spherically- 
averaged auto power spectrum of each CO map using 
Fast-Fourier Transforms (FFTs). We ignore the impact 
of peculiar velocities which have a small impact on the 
spherically averaged power spectra compared to the un- 
certainties in the model. 

Power spectra for several simulated models are 
shown in Figure 2. The figure shows Aq Q co (fc) = 

fc 3 -Pco,Co(&)/(27r 2 ), which denotes the contribution to 
the variance of the CO brightness temperature per ln(fc), 
in units of (/iK) 2 , as a function of k. The strength of 
the CO brightness temperature fluctuations evidently de- 
pends strongly on the model, redshift, and spatial scale. 
At a wavenumber of k = O.lh Mpc -1 and at z = 6.8, the 
amplitude of the power spectrum is A^ co = 0.2(p,K) 2 
for M COjm ; n = 10 8 M Q , while the power spectrum is 
down by about an order of magnitude on this scale by 
z = 9.8. The results appear quite sensitive to M CO:m ; n ; 
in the model with M COjm i n = 10 10 M Q , the power spec- 
trum at k = O.lh Mpc -1 is a factor of ~ 20 smaller 
than in the model with M COim i n = 10 8 M Q . However, 
if we vary M s f im ; n while fixing the SFRD (by adjusting 
the relationship between halo mass and SFR) and fix- 
ing M co min = M s f m i n , this sensitivity is vastly reduced 
(§4-3). ' 

The fluctuations also depend strongly on wavenumber, 
with the fluctuations typically reaching Aqq co (k) ~ 

100(pK) 2 on scales of k <~ a few h Mpc -1 . 

We can understand all of these trends quantitatively 
using a 'halo model' (Cooray & Sheth 2002) type calcu- 
lation. In particular, we expect the auto power spectrum 
of the CO brightness temperature fluctuations to have 
the following form: 



Pco.co(k) = (T co ) 2 



(6) 2 P lin (fc) + 



1 (M 2 ) 



.(14) 



/duty (M) 2 _ 

In this equation Pn n (k) denotes the linear theory density 
power spectrum, while 



<&> = 



A I, 



co. mm 



dMM^b(M) 



M co ,, 



dMM^ 



(15) 



is a mass-weighted bias, which enters here because the 
model CO emission is proportional to host halo mass. 
The second term in Equation 14 involves the mean and 
second-moment of the halo mass function, 



(M 2 ) 
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dMM 2 



co, mm 
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dMM 



Moo,. 



~dM 



(16) 



and represent the shot noise contribution to the power 
spectrum. These equations assume: our usual model 
for connecting CO luminosity and host halo mass, lin- 
ear biasing, that the scale of interest is much larger than 
the virial radius of the relevant halos, and that the halo 
shot noise obeys Poisson statistics. The virial radius of 
a 10 8 M Q halo at z = 7 is 15 co-moving kpc (Barkana 
& Loeb 2001), and so neglecting the impact of the halo 
profile should be a good approximation on the relevant 
scales. The formula also assumes that CO active ha- 
los host precisely one CO luminous galaxy located at 
the halo center - as does our simulated model - likely a 
good assumption for the low mass halos of interest. In 
order to calculate the terms in Equations 14 - 16, we use 
the Eisenstein & Hu (1999) transfer function, which was 
used to generate the initial conditions of the McQuinn et 
al. (2007b) reionization simulations, the Shcth-Tormen 
(2002) mass function, and the Sheth et al. (2001) for- 
mula for halo bias. 

Most of the power spectrum's redshift evolution, and 
its variation with M COjln i n , is driven by its quadratic de- 
pendence on (Too)- (Tco) is hi turn proportional to the 
halo collapse fraction in our model (Equation 13). This is 
partly compensated by (6) 2 , which increases as the host 
halos become rarer and hence more clustered. The de- 
crease in (Tco) is the dominant effect, however, as one 
can see by comparing the M COjm i n = M coo i and 10 10 A/ Q 
curves. The form of Equation 14 also explains the shape 
of the simulated power spectra in Figure 2. On large 
scales the clustering term dominates and the CO fluc- 
tuations trace the underlying density power spectrum, 
while on small scales the Poisson term is the main source 
of fluctuations. The green lines in Figure 2 show that the 
calculations of Equations 14 - 16 agree well with the sim- 
ulation results for the example shown. 11 Note that this 
agreement is expected to the extent that the simulated 
halo mass function and bias match the respective Shcth- 
Tormen (2002) and Sheth et al. (2001) formulas, and 
provided that non-linear biasing and other non-linear ef- 
fects are negligible or swamped by the Poisson term. 

The Poisson term becomes increasingly important as 
the host halos become more massive and less abundant. 
The clustering term also becomes stronger for the more 
massive halos, but the Poisson term grows more rapidly 
with mass, and the two terms become comparable on 
larger scales (smaller k) as the sources become rarer. For 



example, in the M c , 



model at z = 7.3, the 



Poisson and clustering terms are comparable at k slightly 
larger than lh Mpc -1 . On the other hand, if M COjm i n 
is as large as 10 10 M Q , the two terms cross at k « Q.3h 
Mpc -1 (see Figure 2). Note that the models compared in 
the figure have a fixed duty cycle, /duty =0.1. Decreas- 
ing the duty cycle at constant M COjm i n would increase 
the Poisson term at fixed (b) 2 P\i n (k) and Poisson fluctu- 
ations would dominate at still larger scales. Varying the 
CO duty cycle also impacts the mean brightness temper- 
ature, (Tco) ^ /duty, and hence the overall normaliza- 
tion of the model power spectra. The future CO surveys 

11 The slight deficit of simulated power at high k likely results 
because the halos are not a perfect Poisson sample (e.g. Smith et 
al. 2007). 
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Fig. 2. — The auto power spectrum of CO brightness temperature 
fluctuations. The black solid, red dotted, red dashed, blue dashed, 
and cyan dot-dashed curves show simulated CO power spectra at 
different redshifts for various values of M CO In i n . In each case the 
duty cycle is fixed at /duty = 0.1. The green solid line is the 
halo model of Equation 14 for z = 7.3, -M COjIn i n = 1O 1O M , and 
/duty = 0.1. The green dashed line shows the Poisson term, while 
the green dot-dashed curve is the clustering term. 

considered in §7 potentially probe scales in which each 
of the clustering and the Poisson term are important. 

4.3. Parameter Variations 

Before proceeding, let us further examine the impact 
of some of the uncertainties in our modeling. Given the 
success of the halo model of Equation 14 in matching the 
results of numerical simulations, we will use it in this in- 
vestigation. The simplest parameter variation to under- 
stand is the impact of uncertainties in the normalization 
of the Leo — M relation. The average brightness temper- 
ature is proportional to this normalization, and so dialing 
this value up or down simply results in boosting or dimin- 
ishing the strength of the brightness temperature fluctu- 
ations by the normalization squared. Next, increasing or 
decreasing the duty cycle, /duty, also boosts/diminishes 
(Tco) 2 ~~ an d hence the power spectrum normalization 
- as the square of the duty cycle, while simultaneously 
varying the level of Poisson fluctuations as oc l//duty 
For example, for M co , m i n — 1O 8 M the clustering and 
Poisson terms are comparable at k = lh Mpc -1 for 
/duty = 0.1, while these terms are comparable at k = 3h 

MpC" 1 for /duty = 1. 

We have already examined the impact of varying 
-Mco.min in Figure 2; as with the spatially averaged bright- 
ness temperature the power spectrum is less sensitive to 
M sf , min itself provided the SFRD is fixed (§3). This is 
quantified in Figure 3 which shows the impact of increas- 
ing M s f jln i n above the atomic cooling mass, while fixing 



both the SFRD and M r< 



= M si , 



This is accom- 



plished by increasing the normalization of the SFR-M 
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Fig. 3. — The CO auto spectrum for varying M s { m ; n , with 
A^co.min = Af s f jmil> and the SFRD fixed. The SFRD is fixed to 
its value for the model in which the minimum host halo mass is 
the atomic cooling mass. The models are at z = 7.3 and were 
calculated using the halo model. The magenta dashed line shows 
a contrasting model where star formation occurs in halos down to 
the atomic cooling mass, but only halos above M co m j n = 10 10 AfQ 
are CO luminous. 



Fig. 4. — The CO auto spectrum for varying a. The curves fix 
M co m ; n = 10 9 A/q and z = 7.3, while varying the power law index 
a in the CO luminosity-halo mass relationship, Lco °£ M a , at 
constant L(M COim i n ). Boosting a increases the mean emission, the 
clustering strength, and the Poisson fluctuations by giving more 
weight in the CO emissivity to the rarer, more highly clustered 
host halos. 



relation (Equation 5) as M s f jm i n increases. Since (Too) 
varies with the normalization of the SFR-M relation to 
the 3/5th power, this compensates for most of the ex- 
pected drop in (Too)- in addition, raising M s f, m i n in- 
creases both the bias factor and the level of Poisson fluc- 
tuations which further compensates, and actually leads 
to the fluctuations being larger in the high M s f jln i n mod- 
els than in the atomic cooling mass model on some scales. 

Finally we vary the power law index in the Lco ~ M 
relation, Leo oc M a . Our fiducial model adopts a lin- 
ear relationship, a — I, which we choose mainly for 
simplicity. We also assume that the SFR is a linear 
function of halo mass, while z < 3 observations indi- 
cate Leo oc SFR 3 / 5 (§2.2); therefore another sensible 
model is to assume L co oc Af 3 / 5 . Of course, the linear 
SFR-M relation is also chosen for simplicity - regard- 
less, the Leo M 3 / 5 model illustrates the impact of 
weighting the low mass halos more heavily and the mas- 
sive ones less strongly. A contrasting alternative model 
gives the massive halos more weig ht, L co oc M 5 / 3 . It 
is straightforward to generalize the halo model calcula- 
tion of Equations 14-16 for arbitrary power law indices. 
The results of this calculation are shown in Figure 4 for 
z = 7.3, the intermediate case of M COjmin = 1O 9 M , 
with Lqo{M = A/ C o,min) fixed using Equation 12, fixed 
/duty, and each of a = 3/5,1, and 5/3. Increasing a 
increases the contribution of the massive halos, which 
boosts (lco) ! as well as the clustering and Poisson terms. 
The models with a — 3/5 and a = 1 are quite similar, 
and the impact of decreasing a is mostly degenerate here 
with small variations in the Lqq — M normalization. The 



model with a = 5/3, however, produces stronger boosts 
in the average brightness temperature, as well as the bias 
and the Poisson term. In this model, most of the scales 
shown in the figure are dominated by the Poisson term. 
The stronger impact of the a = 5/3 case can be under- 
stood by noting that M l+a dn co /dM oc M~ la+a near 
M = M co . m i n = 10 9 M Q at the redshift of interest. Con- 
sequently, for a = 3/5, 1 the logarithmic contribution to 
the CO emissivity is a monotonically decreasing func- 
tion of mass above M CO;ln j n , but this is not the case for 
a = 5/3. At larger minimum mass and/or at higher red- 
shift, the halo mass function is a more rapidly decreasing 
function of mass, and increasing a has less impact. On 
the other hand, the impact of increasing a is slightly 
larger when decreasing the minimum mass and/or red- 
shift. In the case that the minimum mass is the atomic 
cooling mass, the dependence on a is qualitatively simi- 
lar to the models shown here, but with a more significant 
boost for a — 5/3. 

It would be quite interesting to measure CO auto spec- 
tra of the type shown in Figures 2-4. From sufficiently 
precise measurements one could extract both the cluster- 
ing term and the Poisson term, although this will likely 
be challenging since the power spectrum is fairly close 
to a power law in shape. Nonetheless, it is probably 
reasonable to assume that the power spectrum is dom- 
inated by these terms at large and small scales respec- 
tively. The clustering term is determined by a product 
of the average CO brightness temperature and the bias 
factor of the host halos, while the Poisson term depends 
on the average brightness temperature, and the abun- 
dance/duty cycle of CO luminous galaxies. This would 
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provide valuable information on which types of galax- 
ies are CO luminous at high redshift, how massive their 
host halos are, and what the duty cycle of CO-luminous 
activity is in these halos. Particularly notable is that 
this measurement would constrain the cumulative emis- 
sion from many very faint, unresolved galaxies, which 
individually may be too faint to detect even with deep 
observations from ALMA. More generally, the measure- 
ment would trace the properties of molecular clouds - 
i.e., the locations where stars form - in the very galaxies 
that reionized the Universe. 

It would also be interesting to investigate the angular 
dependence of the CO auto power spectrum. The power 
spectrum will be anisotropic owing to redshift space dis- 
tortions from peculiar velocities (Kaiser 1987), and mea- 
surements of this anisotropy can be used to break model 
degeneracies. A sufficiently precise measurement of the 
quadrupole to monopole ratio of the power spectrum, for 
example, determines the quantity j3 = fl m (z) ' 6 jb « 1/6 
for the CO emitting galaxies. In principle, combining this 
measurement with the spherically averaged power spec- 
trum, allows one to separately determine each of (Too) 
and b from low k measurements, and then infer the Pois- 
son term at high k. 

Yet another interesting quantity to measure is the ratio 
of the power spectrum of fluctuations in the CO (2-1) and 
CO(l-O) lines. The ratio of the root-mean-squared fluc- 
tuations in the two lines gives a (fluctuation-weighted) 
measure of the excitation temperature and probes the 
heating rate in the reionizing sources. 

5. CROSS CORRELATION WITH REDSHIFTED 21 CM 
EMISSION 

Perhaps the most exciting prospect for CO intensity 
mapping, however, is to combine it with future observa- 
tions of the redshifted 21 cm line from the high redshift 
IGM. As motivated in the introduction, the cross spec- 
trum should be less sensitive to systematic effects and 
provide complementary information about the EoR to 
the 21 cm auto spectrum. 

For simplicity, we assume that the spin temperature of 
the 21 cm transition is much larger than the CMB tem- 
perature globally, which is likely a good approximation 
during most of the EoR (Ciardi & Madau 2003, Pritchard 
& Furlanetto 2007). Further, we ignore the impact of pe- 
culiar velocities which should not significantly influence 
our present calculations (Mesinger & Furlanetto 2007). 
With these assumptions, the 21 cm brightness tempera- 
ture at spatial position r can be written as (e.g. Zaldar- 
riaga et al. 2004): 

"1+ l 1 ^ 2 

T 21 {r) = 28mK(xHi)[l + S x (r)][l + S p (r)} -±- 

(17) 

Here (xhi) denotes the volume-averaged neutral hydro- 
gen fraction, S x (r) denotes the fractional fluctuation in 
neutral hydrogen density at spatial position r, and S p (r) 
is the fractional gas density fluctuation. We will also 
use the symbol (x,) = 1 — (%) to denote the volume- 
averaged ionization fraction. The timing and duration 
of reionization are still quite uncertain, and so the red- 
shift at which a given fraction of the IGM volume is ion- 
ized may be different than in our particular reionization 
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Fig. 5. — The cross power spectrum between the CO and 21 
cm brightness temperature fluctuations. Top panel: The absolute 
value of the cross spectrum between 21 cm and CO emission in 
units of (fj,K) 2 at different redshifts and ionization fractions. The 
redshifts at the corresponding ionization fractions are (z, (xi)) = 
(6.90,0.82); (7.32,0.54); (8.34,0.21). The red dashed line adopts 
M m m i n = 10 10 Mq, while the other curves assume that halos down 
to the atomic cooling mass host CO-luminous sources. Bottom 
panel: The cross-correlation coefficient between the two fields as a 
function of wavenumber. 



model. However, Furlanetto et al. (2006b) and McQuinn 
et al. (2007a) show that the size of the ionized regions 
during reionization depend mostly on the ionized frac- 
tion, (xi), rather than the precise redshift at which a 
given volume is ionized. As a result, the shape of the 
cross spectrum at a given ionization fraction is likely a 
more robust prediction than that at a given redshift. 

Using Equation 17, we produce three-dimensional 
maps of the 21 cm field from outputs of the reioniza- 
tion simulations at various redshifts/ionized fractions. 
We then measure the cross power spectrum between 
the 21 cm and CO data cubes as described in the pre- 
vious section. The results of these calculations are 
shown in Figure 5. The top panel shows the abso- 
lute value of the cross spectrum, while the bottom 
panel indicates the cross correlation coefficient between 
the two random fields as a function of wavenumber, 
i"2i,co(k)- The cross correlation coefficient is defined 

1/2 

by r 2 i,co(fc) = f2i,co(fc)/ [Pco,co(k)P2i,2i(k)] and 
is 1 ( — 1) for wavenumbers in which the two fields are 
perfectly correlated (anti-correlated), while it is zero for 
wavenumbers in which the two fields are completely un- 
correlated. 

The simulated cross spectra are similar to those in 
Lidz et al. (2009), and we refer the reader to this pa- 
per for a more detailed discussion, but summarize some 
of the main features here. 12 On large scales, the 21 cm 

12 The main difference with Lidz et al. (2009) is that the CO 



12 



and CO fields are anti-correlated. To understand this, 
consider length scales larger than the size of the ionized 
bubbles during a given stage of reionization. Regions 
that are overdense on large scales contain more galax- 
ies, and are hence brighter in CO emission than typi- 
cal regions. The same regions, however, correspond to 
mostly ionized portions of the 21 cm map, and are con- 
sequently dimmer than average in 21 cm emission. On 
these spatial scales, the two helds are consequently anti- 
correlated. On the other hand, the two fields are uncor- 
related on scales smaller than the ionized bubbles around 
groups of CO-emitting galaxies. This occurs because the 
gas at each point within an ionized region is highly ion- 
ized irrespective of the precise galaxy density. Similarly, 
fully neutral regions do not contain galaxies (unless some 
galaxies have very low ionizing photon escape fractions) . 

The cross-correlation coefficient, T2i,co{k), in Figure 5 
illustrates exactly these trends. On large scales, the sim- 
ulated r 2 i j co(fc) goes to r 2 i,co(fc) = — 1> while it drops 
to zero on small scales. The scale where r2i,co(&) goes 
to zero increases with decreasing redshift as the Universe 
becomes progressively more ionized and the ionized re- 
gions grow. The red dashed line shows that this behav- 
ior is sensitive, however, to the minimum mass of CO- 
luminous galaxies. In the red-dashed line model M co ^ m ; n 
is larger than the hducial value, (M COim ; n = 1O 1O M 
rather than M COjm i n = M coo i), and the cross spectrum 
turns over on larger scales. This happens because the 
more massive halos are more clustered, and hence tend 
to be surrounded by larger bubbles than the less mas- 
sive halos. Although the ionized regions at the stages 
of reionization considered here are much larger than the 
size of HII regions around individual galaxies, it is still 
the case that more massive halos live in larger overdensi- 
ties and tend to be surrounded by larger ionized regions. 
In order to interpret the cross spectrum's turnover scale 
unambiguously, one hence needs to separately constrain 
A^co.minj which may be possible with measurements of 
the CO auto spectrum. 

In principle, the turnover in the cross spectrum be- 
tween a 'traditional' galaxy survey, with resolved galax- 
ies, and the redshifted 21 cm signal, may be easier to 
interpret than the intensity mapping signal considered 
here. However, as discussed in the introduction, tra- 
ditional galaxy surveys are poorly matched to the red- 
shifted 21 cm observations. Furthermore, surveys for 
high redshift LAEs are the only high redshift galaxy sur- 
veys anticipated in the near future with sufficiently wide 
fields of view. The observed abundance of LAEs will be 
modulated by surrounding neutral hydrogen in the IGM, 
which will in turn impact the turnover in the cross spec- 
trum (Lidz ct al. 2009), complicating the interpretation 
of this measurement. Hence the 21 cm- CO cross spec- 
trum may indeed be the most promising approach. 

Let us now consider the amplitude of the cross power 
spectrum, as shown in the top panel of Figure 5. The am- 
plitude of the cross spectrum at k <~ 0.1 hr 1 Mpc varies 
from about 10 2 (^ii^) 2 to about 10 3 (/jK ) 2 , depending on 
redshift and model. These amplitude variations reflect 
several different aspects of the signal, the most impor- 

cmission in our model is proportional to the mass-weighted halo 
abundance, while these authors' galaxy density field is directly pro- 
portional to the halo abundance. 



tant of which are that the cross spectrum is proportional 
to (Tqo), and that the 21 cm signal depends strongly 
on (xi). The dependence on (Tqo) means that the CO 
emission falls off towards high redshift, and with increas- 
ing M COjm i n , although these effects are partly compen- 
sated by increases in (b) with host halo mass and red- 
shift. This behavior is similar to that of the CO auto 
spectrum (§4.2), except that the cross spectrum is linear 
in (Too) and (6), while the auto spectrum is quadratic 
in these quantities. The 21 cm fluctuations are gener- 
ally largest in amplitude on most relevant scales around 
reionization's mid-point, when (xi) — 0.5 (e.g. Lidz et 
al. 2008). This results because the ionized bubbles grow 
larger as reionization proceeds, which initially boosts 
the large scale 21 cm power, while the 21 cm fluctua- 
tions must eventually fall off as the Universe becomes 
completely ionized. In the particular reionization model 
shown, (x^ w 0.5 at z w 7.5, and the cross spectrum 
peaks at this redshift for most relevant wavenumbers. If 
reionization's mid-point occurs at a higher redshift than 
in this model, the peak cross spectrum amplitude may 
be reached at a later stage of reionization - i.e., at a 
larger (x^, (x^) > 0.5. In particular, if reionization's 
mid-point occurs at a sufficiently high redshift, the boost 
in 21 cm power towards (x^ <~ 0.5 may be overcome by 
the fall of in {Too) with increasing redshift. This would 
make it more challenging to detect the cross spectrum 
signal. However, if reionization occurs earlier than in 
our model, the SFRD would also be larger at high red- 
shift and so (Tqo) nray not fall off so rapidly. At any 
rate, this discussion illustrates that measuring the cross 
spectrum amplitude as a function of redshift provides 
additional information, and is another example of how 
the cross spectrum measurement can complement 21 cm 
auto spectrum measurements. 

We now have a sense for the amplitude and shape of 
the 21 cm-CO cross spectrum, and its dependence on 
model, ionized fraction, and redshift. If it is practical 
to measure this signal, it can be used to determine the 
size of the ionized regions around CO-emitting galaxies. 
If the measurement is feasible over a range of redshifts 
and hence at various stages of the reionization process, 
it would provide a direct tracer of the growth of ionized 
regions during the EoR. Combining this with measure- 
ments of the CO auto spectrum can potentially break 
degeneracies owing to uncertainties in the minimum host 
halo mass of CO emitting galaxies, and other uncertain- 
ties in the properties of the CO emitting galaxies. 

6. FOREGROUND EMISSION AND COMBINING 
DIFFERENT CO TRANSITIONS 

Let us move to consider the impact of foreground con- 
tamination, which is one important factor in assessing 
the practicality of measuring the CO auto spectrum and 
the 21 cm-CO cross spectrum. Cleaning the CO fore- 
grounds should proceed essentially as planned for the 
redshifted 21 cm signal: one can distinguish foregrounds 
from the underlying signal by using the fact that the fore- 
grounds are spectrally smooth, while the signal has struc- 
ture in frequency space (e.g. Zaldarriaga et al. 2004, 
Morales & Hewitt 2004, McQuinn et al. 2006, Petrovic 
& Oh 2010). This will separate the CO emission line 
signal from the much larger 'foreground' emission from 
our own galaxy, the primary CMB, point source emis- 
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sion, and other sources. In order to give a rough idea of 
the fidelity at which foregrounds must be removed in the 
CO data relative to the redshifted 21 cm case, it is use- 
ful to consider the ratio of the mean foreground emission 
to the mean signal for each of the 21 cm and CO sig- 
nals/observing frequencies. Adopting 1/iK for the mean 
CO emission signal and 200/j.K for the sky temperature 
at <~ 30 Ghz (e.g. de Oliveira-Costa et al. 2008), the 
ratio of foreground to mean signal is <~ 200 for the CO 
observations, compared to <~ 10 4 for the redshifted 21 
cm observations. From this standpoint, cleaning fore- 
ground emission from the CO observations appears less 
demanding than it does for the redshifted 21 cm signal. 

Note that there is, however, one slight disadvantage of 
cross correlating 21 cm and CO maps as compared to 
the possibility of cross-correlating 21 cm and optically- 
selected galaxies, such as LAEs. The disadvantage is that 
CO measurements at ~ 15 — 30 Ghz will share much of 
the same foreground structure that impacts the 21 cm 
data cubes, unlike the case of the LAEs. For example, 
de Oliveira-Costa et al. (2008) find that a principle com- 
ponent model with only three components can fit maps 
of diffuse Galactic radio emission at frequencies from 10 
Mhz to 100 Ghz, suggesting a highly-correlated origin for 
Galactic foreground emission at these frequencies. 

However, the correlated nature of the CO and 21 cm 
foregrounds is unlikely a big concern. First, cleaning 
the spectrally smooth component of each of the CO and 
21 cm data should remove the contamination extremely 
efficiently. For instance, for the 21cm case, both ana- 
lytic (McQuinn et al 2006) and simulation (Petrovic & 
Oh 2011) calculations show that the foreground residuals 
after cleaning are minuscule, and negligible compared to 
the signal. The main systematic in fact arises from excess 
cleaning, which results in reduction of large scale power. 
A conservative check to eliminate the possibility of corre- 
lated foregrounds can be performed simply by cranking 
up the strength of the foreground cleaning, e.g. by fitting 
the foregrounds with a progressively higher order poly- 
nomial. Moreover, if one can make three-dimensional 
intensity maps in each of the CO(2-l) and the CO(1-0) 
lines for gas emitting at a given redshift, this can provide 
an additional cross-check. For instance, one could take a 
wavelength- weighted difference of maps in the CO (2-1) 
and CO(1-0) lines to remove the highly-correlated fore- 
grounds in these maps, and correlate this difference map 
with the redshifted 21 cm signal. 

Using two rotational lines from gas emitting at a given 
redshift can also protect against the possibility of fore- 
ground 'interlopers' - i.e., line emission from lower red- 
shift galaxies at the observed frequencies of interest. The 
strongest interloper line for the CO(2-l) transition is 
likely to be CO(1-0) emission from gas at lower red- 
shift. Note that contamination from high redshift CO (2- 
1) emission for CO(1-0) intensity mapping during the 
EoR should be negligible. Other possible foreground in- 
terlopers are emission from HCN molecules and radio 
recombination lines, but these are faint compared to the 
CO emission. Regardless, cross correlating data cubes in 
the CO(1-0) and CO(2-l) lines with each other and with 
the redshifted 21 cm data cubes should protect against 
the interloper lines. Visbal & Loeb (2010) consider mul- 
tiple emission line (fine structure lines, CO lines, and 21 
cm) tracers of large scale structure after the EoR and 



consider using cross correlations to evade interloper con- 
tamination in more detail. It would be interesting to in- 
vestigate this further, especially given the uncertainties 
in the ratio of the luminosity in the CO(1-0) and CO(2-l) 
lines, and the potential scatter in this relationship. 

7. DETECTABILITY 

Now that we have theoretical predictions for the CO 
auto spectrum and the 21 cm-CO cross spectrum and 
have discussed systematic effects from foreground con- 
tamination, we turn to consider the detectability of these 
signals. Let us start by considering the CO auto spec- 
trum. First, we consider the CO noise power spectrum 
rather generically, and then we discuss the type of in- 
strument that will be needed to detect the signal more 
concretely (§7.3). 

7.1. Detectability of the CO auto spectrum 

We assume that the CO thermal noise power spectrum 
is a pure white-noise spectrum. The CO brightness tem- 
perature fluctuations are typically of order ~ l(fiK) 2 , 
although with significant dependence on model, redshift, 
and spatial scale. Owing to this, we characterize the 
CO noise by the noise variance in 10 arcminute spatial 
pixels and spectral channels of width Ais/v = 1%; we 
denote this variance by vn and choose values around 
ctat ~ IfiK. The size of the spatial pixels and the spec- 
tral channels are just chosen as convenient numbers at 
which to quote the noise variance; in practice we find 
that slightly smaller pixels are preferable. In the case of 
a pure white-noise random field the noise power spectrum 
may be written as 

-Pn,CO = CjvVpix, (18) 

with V plx denoting the co-moving volume corresponding 
to 1% spectral channels that are 10 arcminutes on a side. 
For reference, this corresponds to 6.75 x 10 3 (Mpc/h) 3 
at z = 7. 

Assuming Gaussian statistics, the variance of the CO 
auto spectrum for a single fc-mode, with line of sight 
component fc|| = \ik and transverse component k\ = 
k ■ ■ isi 

var[P co (k, A1 )] = [Pco(k) + 

P N ,co(fc)e (fc " / ' c "»'" )2+( ^ / ' c — )2 ] 2 . (19) 

The first term on the right hand side of this equation 
is the usual sample variance term, while the second one 
comes from thermal noise in the telescope. The expo- 
nential reflects the limited spectral and spatial resolu- 
tion of the instrument. Here k» res denotes smoothing 
from finite spectral resolution, and k^ res — H(z)/\c{l + 
^K^obs/Az^obs); while k± t . res = 2ir/ (D(z)9 min ) is the spa- 
tial smoothing with D(z) denoting the co-moving dis- 
tance to the redshift of interest, and # m i n giving the an- 
gular size of the spatial pixels. 

In addition, we would like to calculate the variance of 
the spherically-averaged power spectrum. We consider 
logarithmic bins of width e = dink. In this case, the 
minimum variance estimate of the spherically averaged 
power spectrum has a variance of 

1 _ \ " fc 3 Vsurvcy 

a%(k) ^ 4tt 2 a 2 P (k,fi)' [ ' 
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Fig. 6. — Error bar estimates for the CO power spectrum at 
2 = 7. The black squares show error bar estimates for the spheri- 
cally averaged CO power spectrum at z = 7 assuming a theoretical 
model with Af co m i n = ICtMq, while the blue triangles are for 
A^co,min = 10 10 A/q. The CO survey covers a field of view of 25 
deg 2 , with <Tjv = 1^^": spectral channels of width Au = 0.05 Ghz, 
and spans a depth of 68.6 co-moving Mpc/fc (see text). Scales 
roughly to the left of the red dashed line will be lost to foreground 
cleaning, while small scales are lost owing to the limited spatial 
and spectral resolution of the instrument. 



The sum here runs over the upper half-plane (i.e., pos- 
itive n) because we consider the power spectrum of a 
real- valued field, and only half of the Fourier modes are 
independent as a result. It is helpful to note that in 
the case that the noise and signal power spectra are /x- 
independent, the above formula simplifies to Cp(fc) = 
[P(k) + P/v(fc)] 2 /iV Tn , where N m is the number of modes 
in a fc-bin (counting only modes in the upper half plane). 
Here V SUIVcy denotes the co-moving volume of the CO 
survey. If the survey samples a depth AD, centered on 
a redshift z, and covers a field of view on the sky of Jig 
steradians, then Purvey = QsD(z) 2 AD. The sum over fj, 
is restricted by the survey dimensions. 

In order to best measure the 21 cm-CO cross spectrum, 
the CO and 21 cm surveys should aim to have similar 
/c-mode coverage. We anticipate that this requirement 
will fuel the design of the telescopes planned for the CO 
measurement. To achieve this, the CO survey will want 
a wide field of view and several arcminute spatial resolu- 
tion. For our baseline numbers, we assume that the CO 
survey covers 25 deg 2 on the sky, and that AD = 68.6 
Mpc/h, corresponding to a bandwidth of 6 Mhz for the 21 
cm survey we consider shortly. This field of view is com- 
parable to that of LOFAR (Harker et al. 2010), but less 
than the ~ 800 deg 2 planned for the MWA. We assume 
that each spatial pixel has a size of # m i n = 6 arcmin- 
utes, corresponding to k±^ res — 0.58/i Mpc -1 at z = 7. 
For reference, the entire CO survey spans a co-moving 
volume of 2.1 x 10 7 (Mpc//i) 3 co-moving at z = 7. 



Figure 6 shows error bar estimates for two simulated 
models at z — 7. The simulated models have M co , m in = 
10 8 M Q and 10 10 Mq respectively. The power spectrum 
has been averaged in spherical bins of width e = 0.5. As 
mentioned above, the CO intensity mapping experiment 
is assumed to cover 25 deg 2 on the sky, with a thermal 
noise of <jn = lfiK (quoted as the equivalent noise in 10 
arc-minute pixels, and 1% spectral channels), in spectral 
pixels of width Av = 0.05 Ghz, or Av/v = 1.7 x 10~ 3 , 
and 6 arcminute spatial resolution. 13 The results of 
the M CO!in i n = 10 8 Af Q model look encouraging, provided 
these survey parameters are achievable: high significance 
CO power spectrum measurements are possible across 
roughly a decade in spatial scale, k ~ 0.1 — lh Mpc -1 . 
The fractional error bars are larger in the M COim j n = 10 10 
model, which has smaller amplitude fluctuations. This 
quantifies the difficulty in detecting CO fluctuations in 
scenarios in which galaxies with low SFRs are dim in CO. 
If this model is representative, a CO survey with smaller 
thermal noise is required for a firm measurement. Note 
again the distinction between M s { im - m and M COtmin in our 
modeling: if all star-forming galaxies are CO luminous 
- i.e., M CO)m i n = M s f im i n - and the SFRD is fixed, then 
there is relatively little sensitivity to increasing M s f im ; n 
itself (§4.3). If the CO luminosity-halo mass normaliza- 
tion (Equation 12) is a factor of ~ 5 lower than in our 
fiducial model, this would result in a similar reduction in 
S/N as results from increasing M COim ; n from 1O 8 M0 to 
10 10 M Q . 

In the M CO)m i n = 10 8 M© model, the finite spatial and 
spectral resolution prevent measurements on small spa- 
tial scale, while foreground contamination will inevitably 
limit ones ability to measure CO fluctuations on large 
spatial scales. As discussed in §6, one will want to remove 
a spectrally smooth component from the CO data to sep- 
arate the CO emission line signal from the much larger 
spectrally smooth foregrounds. The precise impact of 
the foreground cleaning process will depend on the band- 
width over which the spectrally smooth component is fit, 
the precise cleaning algorithm, and the spectral shape of 
the foregrounds. As a rough guide, it is useful to note 
that if the cleaning is done over a bandwidth correspond- 
ing to a co- moving depth of AD, at least all modes with 
k < 2tt/AD will be strongly impacted by foreground 
cleaning. In fact, modes on smaller scales will have sup- 
pressed power due to aliasing effects, though if correctly 
handled one can still obtain unbiased estimates of the 
power spectrum at these scales, albeit with somewhat 
larger variance (Petrovic & Oh 2011, Liu & Tegmark 
2011). For example, future 21 cm surveys will likely mea- 
sure power spectra over a bandwidth of B ~ 6 Mhz to 
avoid evolution in the signal across the observed band- 
width (e.g. McQuinn et al. 2006). If the foregrounds are 
cleaned over this same bandwidth, this will strongly im- 
pact measurements for modes of length scale larger than 
AD = 68.6 Mpc/ft, or k < 0.092/i Mpc" 1 . This scale is 
indicated by the red dashed line in the figure. One might 
use a larger bandwidth in cleaning the foregrounds, but 

13 For reference, note that if the pixel noise is ctjv = lfJ,K for 
10 arc-minute spatial pixels, and 1% spectral channels, the corre- 
sponding noise variance in 6 arcminutes, Av = 0.05 Ghz pixels is 
ctjv = AfiK. This is obtained assuming pure white noise and using 
Equation 19. 



15 



10 



z 

V} 





1 1 


I 

M . - 10" M c 


CT N = 


0.1 fxK 


«*«.«. = 10'° M s - 








CT N = 


1 fj.K 


X \ 

S \ - 
\ \ 
\ \ 






N \ — 

- — S 

^\ X 








" CT N = 


10 fj.K 











Fig. 7. — Detectability of 21 cm-CO cross spectrum at z = 7, 
(xi) = 0.82 as a function of spectral resolution and pixel noise 
level for the CO measurement. The curves show the total S/N 
(summed over all fc-modes) at which the 21 cm-CO cross spectrum 
may be detected. An S/N of 5 indicates, for example, that the 
cross spectrum can be detected at 5 — a confidence. The solid lines 
are for M cn „ 



10 Mq, while the dashed lines take M co n 



10 iu Af©. 



in this case one will generally require a larger number of 
parameters to adequately describe the foregrounds. This 
will ultimately limit the gain from using a larger band- 
width in the cleaning algorithm (McQuinn et al. 2006). 
Since the CO power spectrum should evolve less with 
rcdshift than the 21 cm power spectrum during reion- 
ization, one can likely extract the CO auto spectrum on 
somewhat larger scales than the 21 cm power spectrum. 
Nonetheless, the dashed red line should provide a conser- 
vative estimate of the scales that are strongly affected by 
foreground cleaning. Measuring the CO auto spectrum 
over even this decade in scale would provide interesting 
constraints on (Too), (b) and the Poisson term, which 
would in turn constrain the star formation rate density 
at high redshift, and the properties of CO emitting gas 
at very early cosmic times. 

7.2. Detectability of the 21 cm-CO cross spectrum 

We now consider the detectability of the 21 cm-CO 
cross spectrum, which is potentially an even more in- 
teresting signal to detect. The variance of the cross- 
spectrum estimate for a single Fourier-mode is a gener- 
alization of Equation 19 (e.g. Furlanetto & Lidz 2007): 

1 2 

var [P 2 i,co(k, n)\ = - [P2i,co(k, fi) + oco(k, ^21 (k, n)\ ■ 

(21) 

In this Equation, o~co(k, (J.) and <7 2 i(fc, (J-) are short-hands 
for the error bar on the CO and 21 cm auto spectra re- 
spectively for a mode with wavenumber fc and /i = fc||/fc, 
while -P2i,co(^j A*) denotes the cross power spectrum. In 



order to evaluate this, we need to first compute the vari- 
ance of the 21 cm auto spectrum which is given by (Mc- 
Quinn et al. 2006): 



var[P 21 (k, (j,)] 



P 2 i(k,M) 



T s 2 ys D 2 AD 
Bti n t n(kx) 



(22) 



The second term in this expression is the 21 cm noise 
power spectrum, while the first term is the usual sample 
variance term. In this equation, T sys is the system tem- 
perature, B is the bandwidth, ij n t is the integration time, 
A is the observed wavelength, A e is the effective area of 
each antenna tile, and n(k±) denotes the number density 
of baselines observing a mode with transverse wavenum- 
ber k± (McQuinn et al. 2006, Bowman et al. 2006, Lidz 
et al. 2007). In order to estimate the error bar on the 
spherically averaged cross spectrum we sum over modes 
in logarithmic fc-bins, as in Equation 20, except restrict- 
ing the survey volume and the sum to modes that are 
common to both the CO and 21 cm surveys. 

We consider the prospects for cross-correlating three- 
dimensional CO maps with 21 cm measurements from 
the MWA. Similar considerations and calculations could 
be performed for the other first generation 21 cm arrays. 
The MWA will consist of 500 antenna tiles each with an 
effective area of A e = 14m 2 at z = 8 (Bowman et al. 
2006). We assume that a fraction of the antennas are 
packed as closely as possible within a 20 m core, and 
that the remaining antennas follow an r~ 2 distribution 
out to a maximum baseline of 1.5 km (Bowman et al. 
2006). We take B = 6 Mhz and t int = 1,000 hrs, and 
T sys = 325 K. The MWA will survey a large field of view 
on the sky of ~ 800 deg 2 ; we assume that only the 25 
deg 2 patch of this covered by the CO survey is available 
for estimating the 21 cm-CO cross spectrum. 

Figure 7 shows an estimate of the total S/N, summed 
over all fc-modes, at which the 21 cm-CO cross power 
spectrum can be detected for various CO survey param- 
eters. These estimates are shown for models with each 
of M CO:inin = 10 8 Af Q and 10 10 M Q . The x-axis shows the 
spectral resolution of the CO measurement, Av/v. For 
a n = 0.1/iif, significant (S/N > 10 — a confidence) cross 
spectrum detections appear feasible in both models. If 
on = 1/iJ^ then a significant detection is only feasible 
in the atomic cooling mass model, while if er/v = 10/iif , 
neither model yields a high S/N measurement. Increas- 
ing the spectral resolution is initially helpful, but the 
gains saturate around Av /v = 0.003 corresponding to 
fell, res ~ 1^ Mpc -1 since the finite spatial resolution 
and thermal noise prohibit measuring higher fc-modes. 
Note that the dependence on spectral resolution is not 
so strong; if the total number of spectral channels that is 
observationally feasible is fixed, it may be advantageous 
to increase the observational bandwidth at fixed spectral 
resolution. 

7.3. CO Instrumental Configuration 

From the above estimates, it appears that significant 
CO auto spectrum detections and 21 cm-CO cross spec- 
trum detections are possible if suitable CO surveys are 
feasible. In the case that a significant fraction of the ion- 
izing sources are CO dim (e.g., M CO:in i n > 10 10 M Q ), a 
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detection will be more challenging, but interesting upper 
limits might be placed. These statements hinge on the 
feasibility of a CO survey covering ~ 25 deg 2 on the sky 
with ~ 6 arcminute, Avjv ~ 0.003 pixels, at a thermal 
noise level of o~n = 0.1 — \\iK (the noise level is quoted 
for coarser 10 arcminute, 1% spectral pixels). 

Here we briefly consider what sort of interferometer 
specifications are needed to match these requirements. 
In order to map 25 deg 2 on the sky in the CO (2-1) line 
at z=7, A bs = 1 cm, antennas of size D ~ 12 cm are 
required. In order to achieve spatial resolution of ~ 6 
arcminutes, maximum baselines of length -D max ~ 6 m 
are necessary. Assuming for starters uniform coverage in 
Fourier space, the thermal noise in an image plane pixel 
is 



1 



°n= .: sy ; / . (23) 

V ZA^Cjnt /cover 

In this equation / COV cr denotes the covering factor of the 
N a antennas, N a D 2 /D 2 aax . We can then estimate the 
covering factor required to achieve cat ~ \\iK in 10 ar- 
cminute, 1% spectral pixels. Assuming T sys = 30 K and 
tint — 1; 000 hrs requires f covcr ~ 1 to reach o~jq = 1//-K"- 
For strictly uniform Fourier coverage this amounts to 
N a ~ 900 antennas, but this is surely a very conservative 
estimate; in reality the antennas will not be uniformly 
distributed and this is likely to boost the sensitivity sig- 
nificantly. Furthermore, the baselines are small enough 
to rotate the entire telescope and help fill in the coverage 
in Fourier space. More detailed investigation is clearly 
warranted, but it seems likely that a few hundred an- 
tenna elements would suffice for the CO measurement. 
An alternative observational strategy is to build a focal 
plane array rather than an interferometer (Bowman et 
al. 2011, in prep.); it would be interesting to compare 
these approaches in detail. 

7.4. 21 cm-CO correlations using next generation 21 
cm Surveys 

Since high redshift CO intensity mapping is a relatively 
new idea (starting with Righi et al. 2008), CO intensity 
mapping instruments may be built on the timescale of 
second generation 21 cm survey instruments. Here we 
briefly consider the prospects for measuring the 21 cm- 
CO cross spectrum using second generation 21 cm sur- 
veys. In order to optimize the S/N of this measurement, 
the two surveys should strive to have similar coverage of 
Fourier modes. These considerations should shape the 
design of each instrument. To name one example, if the 
CO surveys are necessarily limited to fields of view that 
are tens of square degrees, then a 21 cm array design like 
LOFAR which has a smaller field of view than the MWA 
but more collecting area, is advantageous. 

To illustrate the plausible sensitivity of these future 
measurements, we assume a CO survey covering 25 deg 2 , 
with Av = 0.05 Ghz and consider each of o\/v = 0.1/j.K 
and on = lfiK. For the 21 cm survey, we consider an 
array similar to the MWA, except with 5, 000 total anten- 
nas. The antennas are initially arranged as close together 
as possible over an 80 m core, followed by the usual r~ 2 
distribution out to maximum baselines of 1.5 km. This 
is identical to the futuristic version of the MWA dubbed 
'MWA 5000' by McQuinn et al. (2006). We consider the 
ability of these instruments to detect the turnover in the 
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Fig. 8. — Cross correlation coefficient between the 21 cm and 
CO fields with error bars for a futuristic 21 cm survey. The curves 
and error bars are at z = 7. The red dotted line and associated 
1-a error bars assume (x{) = 0.54, while the black solid line shows 
a contrasting model with (xj) = 0.82. The red error bars assume 
that the thermal noise in the CO survey is ajsi = O.l^t-R", while 
the larger blue error bars take crjv = IfiK. The two models can 
be distinguished at high significance. Scales roughly to the left of 
the dashed line are strongly affected by foreground contamination, 
while the error bars blow up on small scales owing to the limited 
spatial and spectral resolution of the CO survey. 



cross-correlation coefficient at z = 7. Here we adopt a 
fiducial model with M s f !m i n = 10 8 M© and (a;,) = 0.54. 
The cross correlation coefficient for this model, with er- 
ror bars is shown by the red dotted line with error bars 
in Figure 8. The black solid line is a contrasting model 
where reionization has progressed further by z = 7 (to 
(x^ = 0.82) and the ionized regions are larger, caus- 
ing the cross correlation coefficient to turn over on larger 
scales in this model. These two models can clearly be dis- 
tinguished at high significance for either level of thermal 
noise in the CO survey, further illustrating the potential 
promise of the 21 cm-CO cross correlation measurements. 
The error bars increase rapidly towards small scales be- 
cause of the limited spatial and spectral resolution of the 
assumed CO survey, and so improvements here could al- 
low even better measurements. 

8. COMPARISON WITH PREVIOUS WORK 

Here we briefly compare with previous work on the 
same topic. The first closely related paper is by Righi 
et al. (2008). Compared to this paper, the main differ- 
ences with the present paper are that we: consider the 
three-dimensional CO power spectrum rather than the 
angular power spectrum as a function of spectral resolu- 
tion, we calculate the 21 cm-CO cross power spectrum, 
and that we consider the detectability of both signals. 
In terms of modeling, Righi et al. (2008) adopt a more 
detailed model for star formation that explicitly ties star 
formation to halo mergers. A second difference, men- 
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tioned already in §2.2, is that they calibrate their as- 
sumed Leo — SFR relation to M82, while we normal- 
ize at significantly lower SFRs. At z = 7, our spatially 
averaged brightness temperature in the CO(2-l) agrees 
with their results to within a factor of a few. Our larger 
brightness temperature is mainly driven by normalizing 
the Leo — SFR relation at lower SFR. The two works 
are in general agreement given the model uncertainties. 

The next related work is Carilli (2011). This work es- 
timated the mean CO brightness temperature by using 
the critical SFRD, along with empirical correlations be- 
tween Leo > Lfir and SFR. Indeed, we followed Carilli 
(2011)'s lead by using the critical SFRD to constrain our 
modeling. Carilli (2011) finds brightness temperatures of 
around 1/iK in the CO(l-O) and CO(2-l) lines (assuming 
the optically thick, high temperature limit). This is in 
good agreement with our own estimates. 14 We expand 
on Carilli's (2011) work by calculating the spatial fluc- 
tuations in the CO emission, the cross power spectrum 
with the redshifted 21 cm signal, and the detectability of 
each signal. 

Finally, the work with the closest overlap to this paper 
is that of Gong ct al. (2011). These authors find a spa- 
tially averaged brightness temperature, in the CO (1-0) 
line, of order \\xK at z <~ 6 — 7. Taken at face value, 
this is broadly consistent with our findings. Interest- 
ingly, however, their assumed Lqo{M) relation is rather 
different than our fiducial model and the models of Righi 
et al. (2008); their model weights the high mass halos 
much more strongly (see their Figure 1). Additionally, 
it appears that /duty ~ 1 in their model in contrast with 
our fiducial model. We have verified that we reproduce 
their mean brightness temperature if we adopt the same 
Leo — M relation. Their Leo — M relation is extracted 
from the Obreschkow et al. (2009) model for CO lumi- 
nous galaxies. We discussed our reservations about using 
this model directly in §2.2. It is not clear what the main 
source of the difference between their Leo — M relation 
and our models are: presumably the Obreschkow et al. 
(2009) model matches the empirical Leo — SFR rela- 
tions at intermediate redshift that are used to calibrate 
our models. The difference then must lie in extrapolat- 
ing this to higher redshift, where semi-analytic models 
have yet to be tested, and to small halo masses that are 
unresolved in the Obreschkow et al. (2009) calculations. 
It would be interesting to investigate this in more detail. 

9. CONCLUSIONS 

We have quantified the possibility of intensity map- 
ping observations in rotational emission lines from CO 
molecules during the EoR. We model both the power 
spectrum of CO emission fluctuations, and the cross- 
correlation between CO and redshifted 21 cm signals. 
These measurements would be extremely interesting for 
three reasons. First, the CO measurement would directly 
inform models for the properties of molecular clouds - 

14 In detail, we make slightly different assumptions in calculating 
the mean brightness temperature and so the near perfect agreement 
with our fiducial model is accidental. Specifically, we conservatively 
assume the same Lqq — SFR relation for the CO(2-l) and CO(1-0) 
lines, but normalize using the sub-linear scaling of Leo with SFR. 
Carilli (2011), on the other hand, assumes a linear Leo — SFR 
relation throughout but uses the optically thick, high temperature 
limit to infer the brightness temperature of each line. Nonetheless, 
the two calculations are in broad agreement given the uncertainties. 



the sites of star formation - in the very galaxies that 
reionize the Universe. Second, the cross spectrum would 
provide a more direct measure of the size of ionized re- 
gions around the ionizing sources than available from the 
redshifted 21 cm signal alone. Third, the cross spectrum 
would help confirm the high redshift nature of a putative 
redshifted 21 cm signal. 

We find that the signal is potentially detectable with 
a relatively modest-sounding interferometer for the CO 
survey. Quantitatively, at z — 7 our fiducial model 
predicts a spatially-averaged brightness temperature of 
{Tco} ~ ^I^K, brightness temperature fluctuations of 
amplitude A co co = 0.2(iiK) 2 at k = O.lh Mpc -1 , 
and a 21 cm- CO cross power spectrum at the level of 
A 2i.co( fc ) ~ 10 3 (a*L:) 2 at k = O.lh Mpc" 1 if the Uni- 
verse is <~ 80% ionized at this redshift. The theoretical 
forecasts are, however, uncertain because we have only a 
limited empirical and theoretical handle on the proper- 
ties of CO emitting gas in the low-luminosity, high red- 
shift galaxies of interest. Here we mention a few possible 
approaches to improving this situation: 

• Early data from ALMA will presumably constrain 
some of the uncertainties in our modeling. 

• Chang et al. (2011) have detected the cross- 
correlation between a post-reionization 21 cm in- 
tensity map and optically selected galaxies from the 
DEEP2 survey at z ~ 1. A powerful proof of princi- 
ple would be to carry out a similar cross-correlation 
between a low-redshift CO intensity map, low red- 
shift 21 cm emission, and optically selected galax- 
ies. In addition to being a strong demonstration of 
the feasibility of these measurements, this measure- 
ment would provide constraints on the aggregate 
properties of CO emission at low redshift where 
star formation is more well understood. It would 
also serve to address any worries about correlated 
foregrounds in the 21 cm and CO maps. A simi- 
lar possibility is to cross correlate low redshift CO 
maps with FIRB maps from Planck and Herschcl 
(e.g. Ade et al. 2011); this would also help con- 
strain which source redshifts dominate the FIRB 
emission. 

• Yet another possible cross-correlation is to combine 
high redshift CO maps with LAE surveys. Since 
the properties of individual LAEs are fairly well 
determined, this will help constrain the nature of 
the CO emitting galaxies in a relatively model in- 
dependent way. 

• Further theoretical work is needed to help under- 
stand plausible values of T ex at high redshift, and 
its dependence on gas density, the CMB temper- 
ature, and on the sources of heat input at high 
redshift. 

It would also be interesting to consider intensity map- 
ping with other emission lines (e.g. Suginohara et al. 
1999, Hernandez-Monteagudo et al. 2007, 2008, Visbal 
& Loeb 2010). The CII 158 fj,m line is the most lumi- 
nous and perhaps most promising other line to consider; 
the advantage of CO is that one can potentially measure 
two rotational lines, providing an important systematics 
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check, and that the necessary millimeter wavelength ob- 
servations for high redshift CII lines are more challenging 
than observations at centimeter-scale wavelengths. 

In conclusion, we believe that CO intensity mapping is 
a potentially powerful approach to measuring large scale 
structure at high redshift. Further theoretical and ob- 
servational work is needed to address the uncertainties 
in predicting the CO signal, but we have outlined sev- 
eral possible approaches to improve this situations. Ulti- 
mately, when used in conjunction with measurements of 
the redshiftcd 21 cm signal, CO intensity mapping may 
provide a powerful probe of the EoR and high redshift 
star formation. 
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